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ABSTRACT 

Within the context of upcoming full-sky lensing surveys, the edge-preserving non- 
linear algorithm Aski (All Sky k Inversion) is presented. Using the framework of 
Maximum A Posteriori inversion, it aims at recovering the optimal full-sky conver- 
gence map from noisy surveys with masks. Aski contributes two steps: (i) CCD im- 
ages of possibly crowded galactic fields are deblurred using automated edge-preserving 
deconvolution; (ii) once the reduced shear is estimated using standard techniques, the 
partially masked convergence map is also inverted via an edge-preserving method. 

The efficiency of the dcblurring of the image is quantified by the relative gain in 
the quality factor of the reduced shear, as estimated by Sextractor. Cross valida- 
tion as a function of the number of stars removed yields an automatic estimate of 
the optimal level of regularization for the deconvolution of the galaxies. It is found 
that when the observed field is crowded, this gain can be quite significant for realis- 
tic ground-based eight-metre class surveys. The most significant improvement occurs 
when both positivity and edge-preserving i\ — £ 2 penalties arc imposed during the 
iterative deconvolution. 

The quality of the convergence inversion is investigated on noisy maps derived 
from the HORIZON-47T N-body simulation with SNR within the range £ cut = 500-2500, 
with and without Galactic cuts, and quantified using one-point statistics (S3 and S4), 
power spectra, cluster counts, peak patches and the skeleton. It is found that (i) the 
reconstruction is able to interpolate and extrapolate within the Galactic cuts/non- 
uniform noise; (ii) its sharpness-preserving penalization avoids strong biasing near the 
clusters of the map (iii) it reconstructs well the shape of the PDF as traced by its 
skewness and kurtosis (iv) the geometry and topology of the reconstructed map is 
close to the initial map as traced by the peak patch distribution and the skeleton's 
differential length (v) the two-points statistics of the recovered map is consistent with 
the corresponding smoothed version of the initial map (vi) the distribution of point 
sources is also consistent with the corresponding smoothing, with a significant im- 
provement when i\ —£2 prior is applied. The contamination of B-modes when realistic 
Galactic cuts are present is also investigated. Leakage mainly occurs on large scales. 
The non-linearities implemented in the model are significant on small scales near the 
peaks in the field. 



1 INTRODUCTION 



In recent years, weak shear measurements have become 
a major source of cosmological information. By measur- 
ing the bending of the rays of light emerging from dis- 
tant galaxies, one can gain some knowledge of the distri- 
bution of matter between the emitter and ourselves, and 
thu s probe the properties and evolu tion history of dark mat- 
ter (|Bartelmann fc Schneiderll200ll ). This technique has led 
to significant results in a broad spectrum of topics, from 



measurements of the proje cted dark m a tter p ower spectrum 
(for the latest results see iFu fc et al. ^OOShh 3D estima - 
tion of the dark matter spectrum ( Kitching et alj |2006| ). 
studies of the higher order moments of the dark matter 
distributio n, selection of source candidates for subsequent 
follow-ups (jSchirmer et alj | 2007f). and recons truction of the 
mass distribution fro m small (|jee et alj|2007h to large scales 
jMassev et all 120071 '). In view of these successes, numerous 
surveys have been planned specifically to use this probe 
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either from ground-based facilities (eg VST-KID^], DES^ 
Pan-STARRffl LSST0) or space-based observatories (EU- 
CLIE0, SNAffl and JDE^ll)- More generally, it is clear 
that weak lensing will be a major player in the future, as 
it has been identified by different European and US working 
groups as one of the most efficient way of studying the prop- 
erties of dark energjjf). Data processing is an important issue 
in the exploitation of weak lensing of distant galaxies. The 
signal comes from the excess alignment of the ellipticities of 
the observed galaxies. Assuming one can ig nore or deal with 
spuri o us alignments due t o intrinsic effects ftlirata fc Seliakl 
12004 lAubert et all |2004| ; iPichon fc Bernardeau jl999l 'l. or 
due to spurious lensing effects jBridle fc Abdallall2007r i. the 
weak lensing signal will thus come from a small statistically 
coherent ellipticity on top of the random one of each object. 
Any result obtained with weak lensing on distant galaxies 
is thus conditioned by the quality with which shape param- 
eters of the galaxies are recovered. This issue has of course 
been raised by the weak lensing community and tackled by 
the SHear Test ing Program working group (|Massevl [20071 ; 
iHevmansI [20061 ') whose effort have allowed for a fair com- 
parison of the existing techniques. Schematically, the mea- 
surement of the shape parameters of the galaxies can be 
seen as a two-step process. First, one must correct for the 
non-idealities of the images due to atmospherical seeing (for 
ground-based telescopes), and telescope and camera aber- 
rations. Indeed, these effects translate into an asymmetrical 
beam, which is varying between two images, and even pos- 
sibly in the field of one image. Typically, the asymmetry 
induced by the instrumental response is much larger than 
the ellipticity to be measured. After this preprocessing step, 
a shape determination algorithm can be applied, and some 
estimation of the ellipticity of the object recovered. Stars, 
defects in the images, and objects too close to each other 
after deconvolution have to be removed from the final cat- 
alogue so as to avoid contamination from erroneous shape 
measurements. 

After these operations, one obtains a catalogue of posi- 
tion and shape parameters. Many techniques exist for recov- 
ering the weak shear signal from this catalogue. For exam- 
ple, a lot of efforts have been devoted to the measurement of 
the shear two-point functions. The most used method is the 
two-point functions; however measurement of the so called 
Mass Aperture averaged two-point function, which is the re- 
sult of the convoluti on of the shear two-po int functions by 
a compensated filt er (ISchneider et al.|[2002l ) is becoming the 
preferred method (|Fu fc et al. 20081 ) . This scheme includes 
the separation between the curl-free convergence-field two- 
point function, and the residual curl mode that can arise 
from incomplete PSF correc tion or intrinsic galaxy align- 
ment l|Crittenden et al1l2002l ). For three-point functions, dif- 



1 http : //www . astro-wise . org/projects/KIDS/ 

2 https : / /www . darkenergy survey . org/ 

3 http : //pan-starrs . if a.hawaii . edu/ 

4 http://www.lsst.org/ 

5 http://www.dune-mission.net/ 

6 http://snap.lbl.gov/ 

7 http : //universe .nasa. gov/program/probes/jdem. html 

8 see, on the European side http : //www. stecf . org/coordination/ 
and on the US side, http://www.nsf.gov/mps/ast/aaac.jsp 
and http: //www .nsf .gov/mps/ast/detf . asp 



ferent resummation schemes have b een proposed, either us- 
ing direct measurement of the s hear l|Bernardeau et al]|2002l ; 
iBenabed fc Scoccimarro] 120051) or using the Mass Aper ture 
filter (|Takada fc Jainll2003l ; iKilbinger fc Schneider!l2005h . 

Other applications (source detection and fit, some to- 
mography algorithms) call for an estimation of the map of 
the convergence field. A convergence map can also be used to 
measure the two- and three-point functions as well, even if, 
as we will see later this is not optimal. For these reasons an 
important amount of work has already been devoted to the 
reconstruction of the convergence map (Ivan Waerbeke et alj 
Il999l ; ISeitz et all 1 19981 ; iBartelmann et al.lll996t ). The prob- 
lem in this reconstruction lies in the inversion of the non- 
local equations linking the convergence field «, and the ellip- 
ticities of the galaxies, while controlling the noise and avoid- 
ing pollution from the spurious curl modes. Moreover, even 
assuming that the ellipticity catalogue was a noise free es- 
timation of a curl-free underlying shear, the inversion could 
only be exact up to a global translation given the functional 
form of the equation. Thus Bayesian techniques that use 
a priori properties on the solution to regularize the inver- 
sion problem are well suited to the reconstruction of k. Pre- 
vious works on the topic have explor ed different sets of a 
priori and regularization techniques dMarshall et al. 1200 2| ; 
IStarck et al1l2005l ; ISeitz et al.lll998l ; iBridle et al.lll99sh . The 
primary goal of those investigations being the measurement 
of the mass distribution in clusters, most of them are dealing 
only with finite regions of the sky. For the same reason those 
papers have been extended to include strong-lensing effects 
that can be observed around the cluste r whose mass is being 
reconstructed using their lensing effect ([Cacciato et alj2006l : 



iBradac et all 120051 ; lHalkola et all l200rj ; |jee et al.ll2007l ) 

In this paper, we will focus on the optimal reconstruc- 
tion of the k Held from very large, and possibly full-sky maps, 
of the sky. We will thus only be interested in the weak lensing 
regime including the onset of the quasi-linear regime, where 
the non-linearities of the relation linking the ellipticities of 
the galaxies to the shear cannot be safely neglected. We 
will propose a self calibrated regularization technique, that 
can be co mpared to multi resolution methods or wavelet 
approach |Starck et al.ll2005l ; lAbrial et al.ll2008r i. and use a 
ii — £2 regularization scheme to perform a sharp feature 
preserving inversion. One of the biggest issues we will have 
to cope with is the incomplete coverage of the sky. We will 
show how our technique can deal with irregular coverage and 
masked portions of the sky. 

Specifically, Section [2] shows how self calibrated non- 
parametric £\ — £2 deblurring can improve the construction 
of reduced shear, hence convergence maps. Section [3] de- 
scribes the model for the reduced shear, the corresponding 
inverse problem, and the optimization procedure. Section 2] 
investigates the quality of the global reconstruction; in par- 
ticular, it probes the asymmetry /kurtosis of the recovered 
maps, its topology (total length and differential length of 
the skeleton), the recovered power spectra, the point source 
catalogue with and without galactic star cut. The leaking 
of B-modes induced by the Galactic cut is also investigated. 
Finally, Section [5] discusses implications for upcoming full- 
sky surveys and wraps up. 

Appendix [X] describes the star removal algorithm (imple- 
mented for the cross validation estimation of the optimal 
level of smoothing required to deconvolve the crowded im- 
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ages), Appendix [B] details the k inverse problem on the 
sphere while Appendix[C]derives the local plane correspond- 
ing approximation. Appendix [D] describes the construction 
of realistic k maps from large N-body simulations. 



2 DEBLURRING OF CROWDED FIELDS 

The first step involved in reconstructing a full-sky map of 
the convergence on the sky requires estimating ellipticity 
and orientation maps from wide angle CCD images of large 
patches of the sky. Whether the experiment is ground-based, 
or space-born, it is advisable to correct for the effect of the 
instrumental response, in particular when mapping more 
crowded regions closer to the galactic plane. Indeed, the 
PSF-induced partial overlapping of galaxies within the field 
of view will bias the estimation of the reduced shear. What 
we will describe here would correspond to a method belong- 
ing to the "or ange" q u adran t of the classification proposed 
m table 3 of |m assc v (2007). Current methods have been 
designed for deblurring of isolated objects and are conse- 
quently less efficient in deblurrinng blended objects. As 
a first step towards building a full-sky map maker, let us 
therefore address the issue of deblurring crowded fields via 
regularized non parametric model fitting, and assess its ef- 
ficiency in the weak lensing context. 

In particular we will show that cross validation as a 
function of the number of stars removed yields an automatic 
estimate of the optimal level of regularization for the decon- 
volution of the galaxies. When the observed field is crowded, 
this gain can be quite significant for realistic ground-based 
eight-metre class surveys. The most significant improvement 
occurs when both positivity and edge-preserving l\ — £2 
penalties are imposed during th e iterative deconvolution. 

2.1 Deblurring as an inverse problem 

2.1.1 Regularized solution 

Since observed objects are incoherent sources, the observed 
image depends linearly on the sky brightness distribution: 

y(uj) — J /i(w, w') x(w') d«' + e(«) , 

where y(u) is the observed distribution in the direction uj, 
h(uj,u}') is the atmospheric and instrumental point spread 
function (PSF) which is the distribution of observed light in 
the direction lj due to light coming from direction u' , x(u>') 
is the true sky brightness distribution and e(uj) is the noise. 
After discretization: 

y = Ux + e, (1) 

where y is the vector of pixel intensities in the observed 
image (the data), H is the matrix which accounts for 
the atmospheric and instrumental blurring, x is the (dis- 
cretized or projected onto a basis of functions) object 
brightness distribution and e accounts for the errors (pixel- 
wise noise and modelisation approximations) . Deblurring re- 
quires estimating the best sky brightness distribution given 
the data. Since the atmospheric and instrumental PSF re- 
sults in a smoother distribution than the true one, it is 
well known that de-blurring is an ill-conditioned problem 
( (|Richardsonlll972l : ISkilling et ai1ll979l ; iTarantola fc Valettei 



1 19821 : I Pichon fc Thiebautl ll99Sl ; IPichon et all l200ll )). In 

other words, straightforward deblurring by applying H 1 
to the data y would result in uncontrolled amplification of 
noise: a small change in the input data would yield unac- 
ceptably large artifacts in the solution. Regularization must 
be used to overcome ill-conditioning of this inverse problem. 
This is achieved by using additional prior constraints such 
as requiring that the solution be as smooth as possible, while 
being still in statistical agreement with the data and while 
imposing that the brightness distribution is positive. Fol- 
lowing this prescription, the Maximum A Posteriori (MAP) 
solution x M is the one which minimizes an objective function 
Q(x): 

= argmin Q(x) , with: Q(x) = C(x) + fi1Z(x) , (2) 

where £(x) is a likelihood penalty which enforces agreement 
of the model with the data, TL(x) is a regularization penalty 
which enforces prior constraints set on the model, and /i > 
is a so-called hyper-parameter which allow the tuning of the 
relative weight of the prior with respect to the data. Hence 
the MAP solution is a compromise between what can be 
inferred from the data alone and prior knowledge about the 
parameters of interest. Assuming Gaussian statistics for the 
errors e in equation |T}, the likelihood penalty writes: 

£(x) = (H • x - y) T .W-(U-x-y), (3) 

where the weighting matrix W is equal to the inverse of the 
covariance matrix of the errors: W = Cov(e) -1 . 

The most effective regularization for ill-conditioned 
problems such as deconvolution of b lurred images c onsists 
in imposing a smoothness constraint (jThicbaut 2 0051 ) . Then 
the regularization penalty writes: 

n(x) = Y i <f(Ax j ), (4) 

i 

where Axj is the local gradient of x and (j> is some cost func- 
tion. The local gradient of x can be approximated by finite 
differences: Ax = D ■ x where D is a linear finite difference 
operator. For instance, in 1-D: Axj = (D • x)j = Xj+i — Xj. 
To enforce smoothness, the cost function (j> must be an in- 
creasing function of the magnitude of its argument. Very 
common choices for <j> are: the £2 norm, the l\ norm, or an 
i\ — I2 norm. For our deblurring problem, we have consid- 
ered different priors (quadratic or l\ — I2 smoothness) pos- 
sibly with an additional positivity c onstraint. We have used 
generalized cross validation (GCV. (|Wahbal |l99Ch ) applied 
to the circulant approximation of the quadratic problem to 
estimate the optimal regularization level fi. These different 
possibilities and their effects on the recovered images are 
discussed in details in what follows. 

Finally, to solve for the constrained optimization prob- 
lem we used the v mlmb algorithm from Optim- 
Pack ((thiebautI I2002T ). Vmlmb (for Variable Met- 
ric, Limited Memory, B ounded) makes use of a BFGS 
jNocedal fc Wright! [2006 ) update of the approximation of 
the Hessian (matrix of second partial derivatives) of Q(x) to 
derive a step to improve the parameters at every iteration. 
This strategy only requires computing the objective func- 
tion, Q(x), and its gradient (partial derivatives) V x Q(x) 
with respect to the parameters x. The BFGS update is lim- 
ited to a few last steps so that the memory requirements 
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Hyperparameter by sub-pixel CLEAN + GCV 



2.1.2 Quadratic regularization and Wiener proxy 




number of removed point-like objects 

Figure 1. Hyper-parameter chosen by GCV as a function of 
the number of stars removed by our star removal algorithm (Ap- 
pendix fXJ. Note that this curve reaches a maximum correspond- 
ing to the moment when all the stars have been removed. Indeed 
stars correspond to high frequency correlated signal, while the 
wings of galaxies (for which the core has been erroneously re- 
moved) also give rise to such signal. In between, when all stars 
have been removed, while no galaxies has yet been deprived of 
its core, the amount of correlated high frequency signal reaches a 
minimum, or equivalcntly the GCV estimated value of fi reaches 
a maximum. 



^-sidc 


128 


256 


512 


1024 


2048 


time for one step (s) 


0.13 


0.59 


2.13 


8.48 


34.3 


number of steps (s) 


13 


12 


9 


13 


24 


total time (s) 


2.6 


10.4 


33.4 


171.1 


1129.3 



Table 1. the performance of the optimization of the linearized 
inversion problem ra s jd c >S'Fs lt as a function of n s id c for an octo 
Opteron in OpenMP. 



remains modest, that is a few times the number of sought 
parameters, and the algorithm can be applied to solve very 
large problems (in our case, there are as many parameters as 
the number of pixels in the sought image). Finally, Vmlmb 
accou nts for bound constraints by means of gradient projec- 
tions l|Nocedal fc Wrightj|2006r i. For a convex penalty Q(x), 
VMLMBis guaranteed to converge to the unique feasible min- 
imum of Q(x) which satisfies the bound constraints; for a 
non-convex penalty, VMLMBbeing based on a descent strat- 
egy, it will find a local minimum depending on the initial set 
of parameters. 



Using the finite difference operator D and an I2 norm for 
the regularization and ignoring for the moment the posi- 
tivity constraint, the MAP solution is the minimum of a 
quadratic penalty which simply involves solving a (huge) 
linear problem: 

cc M = argmin {(H • x — y) ■ W • (H • x — y) 

+li(D-xf ■ (D-aj)j 

= (H T ■ W H + AtD 1 ' • d) H t W y, (5) 

providing the Hessian matrix H T • W -H + /i D T • D is non- 
singular, which is generally the case for [i > 0. Owing to the 
large size of the matrices involved in this equation (there 
are as many unknown as the number of pixels), the linear 
problem has to be iteratively solved (by a limited memory 
algorithm such as vmlm) unless it can be diagonalized as 
explained below. The solution, equation ([5}, involves at least 
one parameter, /x, which needs to be set to the correct level of 
regularization: too low would give a solution plagued by lots 
of artifacts due to noise amplification, too high would result 
in an oversmoothed solution with small details blurred. The 
optimal level of smoothing can be computed by generalized 
cross vali dation (GCV) by minimizing wit h respect to /x the 
function jGolub et al.lll979l : IWahballl990h : 



GCV(^) = 



(Ag ■ y - yf ■ W ■ (Ag ■ y - y) 

[1 - tr(A M )/Af 



where N is the number of data (size of y ) and A M 
Xfj) is the so-called influence matrix, in our case: 



Au = H • f H • W ■ H 



(iD T ■ D 



(6) 
V„(H- 

■ (7) 



Computing the value of GCV(p) involves: (i) solving the 
problem to find the regularized solution x M and compute 
A M y — H • cc M ; (ii) estima te the trace o f A M perhaps by us- 
ing Monte Carlo methods (|Girardlll989l ) since the influence 
matrix is very large. The computational cost of stages (i) 
and (ii) is similar to that of a few solvings of the quadratic 
problem. Since this has to be repeated for every different 
value of the regularization level, finding the optimal value 
of /x by means of GCV can be very time consuming unless 
the problem can be approximated by a diagonal quadratic 
problem (for which matrix inversions are both fast and triv- 
ial). 

For this purpose, we introduce the proxy problem cor- 
responding to white noise and circulant approximations of 
the operators H (convolution by the PSF) and D (finite 
differences). Then the weighting matrix becomes: 



Wi, 



Si , 3 /a 



where a = Var (m). 



where a 2 = Var(e^) is the variance of the noise. In the special 
case where the PSF is shift-invariant, H is a convolution 
operator which can be approximated by a block Toeplitz 
with Toeplitz block matrix that can be coputed very quickly 
by means of FFT's: 



Ha; 



diag(F • h) ■ (F ■ x) 



(8) 



where h is the PSF (the first row of H), F is the forward 
DFT operator, and diag(v) is the diagonal matrix with 
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its diagonal given by the vector v. This discrete convolu- 
tion equation assumes that F u j = exp(— 2 i ir u n j n /N„) 
where N„. is the length of the n th dimension, j„ = 
0, . . . , N n — 1 and u n — 0, . . . , N n — 1 are the indices of the po- 
sition and discrete Fourier frequency along this dimension. 
In this case, the inverse DFT is simply F" 1 = F H /N to t with 
iVtot the total number of elements in x and the H exponent 
standing for the conjugate transpose. With these approx- 
imations and definitions of the DFT, the likelihood term 
writes: 

£{ x ) = — ll H • x - yf - M 1 2 J2 ft™ - M 2 > ( 9 ) 

u 

where h u is the transfer function (the DFT of the point 
spread function) and y u and x u respectively the DFT of the 
data and of the sought image. Note that the exact normaliza- 
tion factor, here l/JV to t, depends on the particular definition 
of the DFT. 

Similarity, ignoring edges effects, the finite difference 
operator D along n th direction can be approximated by: 



with 



D, 



diag(d„) • (F • x) 



where d n is the DFT of the first row of D„ 
quadratic regularization writes: 



(10) 
then the 



TZ(x) 



with: 



ID-asI 



Ni 



- ]Tr u \x u \ 2 , (11) 



(12) 



for first order finite differences and our choice for the DFT. 
Note that any r„ ^ being an increasing function of the 
length |u| of the spatial frequency could be used instead and 
would result in imposing a smoothness constraint although 
with a different behaviour. Finally putting all these circulant 
approximations together, the quadratic problem to solve is 
diagonalized in the DFT space and trivially solved to gives 
the DFT of the MAP solution: 



K Vn 



\h u \ 2 + fj,a 2 r v 



(13) 



the asterisk exponent denoting the complex conjugate. Note 
that this circulant approximation of the solution is very fast 
to compute as it involves just a few FFT's. This expres- 
sion of the MAP solution is very similar to what would give 
Wiener filter which would be exactly achieved by setting 
the term fi r u equals to the reciprocal of the expected image 
powerspectrum in equation (|13[) . Since, in our case, the im- 
age powerspectrum is unknown a priori, we have to choose 
the regularization shape r u and derive the optimal level of 
smoothing, for instance, by means of GCV. Thanks to the 
circulant approximation made here, GCV criterion is now 
very easy to compute as: 



~ F 1 • diag(a M ) ■ F , with: a M ,« 



|/t u | 2 + fia 2 r u 



and tr(A M ) = J2 U <ty,u/Wtot, hence: 

Mot E^jU l£"l 2 



GCV(/i) = 



V 2 E„ *M,«] 



(14) 



= 1 



MO" r u 



+ jj,a 2 r n 



(15) 



In practice, for the optimization of equation equa- 
tion (1131) is taken as a starting point together with the choice 
of fi given by the minimum of equation (|14[) . Then the opti- 
mization of equation ([2} is carried with possibly non station- 
ary weights, while iterating back and forth between model 
and data space. 



2.1.3 Crowded fields and star removal 

Even though the estimation of the ellipticities does not re- 
quire per se the deconvolution of the galaxies, it is shown 
below that this estimation is significantly improved by de- 
convolution when the fields of view are crowded and polluted 
by foreground stars: indeed galaxies and stars overlap less 
when deconvolved, which reduces the fraction of erroneous 
measurements. Unfortunately, when these stars are present, 
they significantly bias the estimation of the hyper parame- 
ter, fj,, since stars correspond to high frequency correlated 
signal which leads to an underestimation of the optimal level 
of smoothing (for the galaxies) by cross validation. This is 
best seen in Figure [l] which displays the evolution of the 
hyper-parameter which minimizes GCV as a function of the 
number of stars removed by our star removal algorithm, see 
Appendix [A] Interestingly, it suggests that GCV could be 
used as a classifier. 



2.1.4 £1 — £2 penalty and positivity 

The drawback of using a quadratic (£2 ) norm in the regular- 
ization is that it tends to over-smooth the regularized map 
especially around sharp features as point-like sources (i.e. 
stars) and the core of galaxies. This is because the regular- 
ization prevents large intensity differences between neigh- 
boring pixels and result in damped oscillations (Gibbs ef- 
fect). Such ripples hide any faint details in the vicinity of 
sharp structures. To avoid this, it would be better to use a 
regularization which smoothes out small local fluctuations 
of the sought distribution (here the deblurred image), pre- 
sumably due to noise, but let larger local fluc tuations arise 
occasionally fsee lAubert fc Kornprobstl (|2008l ) and reference 
therein). This can be achieved by using a &\ — £2 cost func- 
tio n <j) in equation Q. A possible £\ — £2 sparse cost function 
is (|Mugnier et alj|2004 ): 



2e' 



log 



1+ \ r -\ 



For a small, respectively large, pixel differences r, 
the following behavior 



(16) 
(r) has 



<t>{r) 



r when |r| <C e . 

2 \er\ when \r\ >• e , 



which shows that, as required, the £\ — £2 penalty behave 
quadratically for small residuals r's (in magnitude and w.r.t. 
e) and only linearly for large r's. The derivative, needed for 
the optimization algorithm, of the £1 — £2 penalty writes: 

2er 



0'M = 



e + r 
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Figure 2. An example of virtual fields generated with SKYMAKER to be fed to SEXTRACTOR before and after deconvolution using the 
different rcgularizations described in the text. From top to bottom and left to right, a galaxy field image and the corresponding "true" 
field, a galaxies field with stars and a crowded galaxy field with stars (10 6 stars/arcmin 2 ). The exposure time is fO seconds and the 
seeing is 1" for the VLT with VfMOS. The background field corresponds to the actual size of the corresponding observed images. 



An additional possibility to improve the restitution of faint 
details with level close to that of the background is to ap- 
ply a strict positivity constraint. This is achieved by using 
vmlmb, a modifie d limited memory variable metric method 
(|ThiebautI|2002| ). which imposes simple bound constraints 
by means of gradient projection. This yields a reduction of 
aliasing by bounding the allowed region of parameter space 
which can be explored during the optimization. 



2.2 Numerical experiments 



The public package SkyMaker (|Erben et AL.IIioOlT ) was 
used to generate galactic and stellar fields from ellipticity 
and magnitude catalogues. Table [2] summarizes the main 
parameter corresponding to the VLT with a VIMOS instru- 
ment, a worse case situation compared to upcoming space 



missions. 



A regular grid of 12 x 12 galaxies of magnitude 20 with 
random orientation is produced twice (with the same ran- 
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Figure 3. the error in ellipticity as a function of the elliptic- 
ity (measured by SEXTRACTOR) for a set of 50 images (such as 
those shown on Figure [2]l cither directly on the image (medium 
diamonds), deconvolved with £2 gradient penalty function with 
enforced positivity (light squares) and £\ — £2 gradient penalty 
function with positivity (dark circles). For each set, the elliptic- 
ity is also measured directly on the raw image. Note that, as ex- 
pected, the error on the bias is largest for circular galaxies, since 
deconvolution will tend to over amplify departure from circular 
symmetry. 



dom seed), one corresponding to a fixed seeing and a given 
exposure time, while the other assumes zero noise and zero 
seeing for a set of 512 x 512 pixels images, see Figure [2] 

The background level and the amplitude of the back- 
ground noise is first estimated automatically from the 
histogram of the pixel va lues and fed to Sextractor 
(|Bertin fc ARNOUTslll99ot ) which then estimates the po- 
sition, the flux, the orientation and the ellipticity for all the 
galaxies in the field. Here the ellipticity is defined as 1 — b/a, 
where a and b are the long and short axis. This procedure 
is reproduced 50 times with different realizations. The mea- 
sured and the recovered ellipticity are compared, together 
with flux and orientation for all the galaxies in the field. In 
this set of simulations the prior knowledge of the position 
of the galaxy is used to minimize errors which might arise 
while using SEXTRACTOR: the recovered galaxy is chosen to 
be that which is closest to the known input position. The 
median and interquartile of the error (difference between the 
"true" and recovered) in ellipticity versus the ellipticity is 
computed for a range of exposure time; this procedure is 
iterated for the three deconvolution techniques used in this 
paper (Wiener, £2 with positivity, £\ — I2 with positivity). 
An example of such a plot is shown in Figure [3] Clearly the 
bias in the recovered ellipticity increases with the ellipticity 
and the amount of noise in the image (via poorer seeing or 
shorter exposure time). As expected, the Wiener deconvo- 
lution is the least efficient of the three methods, since the 



linear penalty does not avoid some level of Gibbs ringing. In 
contrast the £2 penalty with positivity avoids partially such 
ringing, while the l\— 12 penalty works best at recovering the 
input eccentricity with a consistent level of bias below 10 % 
for an ellipticity in the range [0.1, 0.8[. Note that this bias is 
relative, not absolute. If an alternative shear estimator that 
doesnt consider deconvolution is accurate to a level of, say 
1%, the expected bias after deconvolution will be below 0.1 
%. 

Interestingly, there is also a residual bias (even for 
longer exposure times) for small ellipticity galaxies, which 
arises because noise induced departure from sphericity is 
amplified by the deconvolution. Note that the Wiener de- 
convolution is significantly faster than the iterative decon- 
volution with positivity (with £2 or £\ — £2 penalties). Posi- 
tivity improves significantly the deconvolution, but will de- 
pend critically on the ability to estimate the background. In 
the present simulations, the level of background is automat- 
ically estimated while looking at the histogram of the pixels. 
Finally the £\ — £2 regularization significantly improves the 
restoration of fields of stars and galaxies, because the stars 
and the cores of galaxies are very sharp. These non- linear 
iterative methods are slower than the Wiener filtering, but 
can account at no extra cost for non uniform noise, or satu- 
ration and masking. Their convergence can be considerably 
boosted when they are initiated by the Wiener solution. 

For any such plot, two numbers are defined which sum- 
marize the trend. The mean error (averaged over the various 
ellipticities) e, and the mean of the interquartile, Ae were 
measured. The quality factor, QF is defined to be the ratio 
of the sum of this mean error and the mean interquartile 
for the image without deblurring, divided by the sum of the 
mean error and the mean interquartile for the deconvolved 
image for the three techniques (Wiener, £2 and £1 — ^2). This 
reads 



^method ~t~ A6 me ^l 1 od 



The evolution the quality of the ellipticity measurement is 
traced versus seeing conditions and signal to noise (exposure 
time) in two regimes: a galaxy-only field, and a galactic field 
with a crowded star content where the number of stars per 
square degree reaches 10 s stars/arcmin 2 . These two regime 
represent high and low Galactic region respectively. Figure 
[4] displays the evolution of QFWiener (diamonds), QFe 2 (tri- 
angles), and QFt 1 _e 2 (circles), as a function the exposure 
time of 1, 10, 100 and 1000 seconds respectively, and two 
seeing conditions of 1.2" and 0.7". No stars are present in 
the field on the left panel of Figure 3J whereas its right panel 
displays the three QF estimators for a field with a realistic 
10 stars per square degree. Aski achieves efficient deblur- 
ing in this regime. It remains to be shown that regularized 
deconvolution obtained through (sparse) parametric local 
decomposition of both PSF and objects (as done e.g. with 
shapelet-based methods) can properly deblur blended ob- 
jects. 

Now that we have shown that state of the art automated 
positive edge-preserving deconvolution of deep sky images is 
mandatory to get good quality shear estimates (most impor- 
tantly in the context of crowded fields) , let us conclude this 
section by a leap forward, and assume from now on that we 
have access not only to discrete measurements of elliptici- 
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Figure 4. Left panel: the relative error quality factor (see main text) as a function of the log exposure time for the three methods, 
respectively Wiener filtering (diamonds), £2 gradient penalty function with enforced positivity (triangles) and £\ — £2 gradient penalty 
function with positivity (circles). Two seeing conditions are investigated, corresponding to a good (0.7") and a fair (1.2") seeing condition. 
These simulations assume that no star are present in the field, and correspond to a set of non overlapping galactic disks with random 
orientation and magnitude 20 in V (see Fig- - T ne telescope setting correspond to the VIMOS instrument on an 8 meter VLT. Right 
panel: the quality factor as a function of the log exposure time, but this time while allowing for stars in the field. The star count is 10 5 
stars per arcmin 2 . As discussed in the text, the penalty weight is estimated via generalized cross validation (GCV) on a temporary image 
where all stars are automatically removed via blind cleaning as describedin Appendix A. Here the removal of stars is essential since the 
GCV hyper-parameter (which sets the level of smoothing in the deconvolved image) varies by orders of magnitudes in the process (see 
Fig [Tj for a discussion) and would be otherwise underestimated. 



ties over a significant fraction of the sky, but also that this 
point like process has been re-sampled. Indeed, since it is 
beyond the scope of this paper to carry out a full-sky decon- 
volution and reconstruction at the resolution of 0.7" (This 
would amount to about 10 12 pixels!), it is assumed from now 
on that a full-sky catalogue of vector reduced shear exists 
and that the interpolation/re-sampling of the correspond- 
ing map on a uniform grid over the sphere has been done, 
together with an estimate of the corresponding shot noise. 
In other words, we skip the critical step of optimal shear 
estimation, which has already been addressed by the STEP 
man sll200fj ; |Mass"evll2007l ) working group. In this paper, 
we extract the virtual catalogue from a state of the art simu- 
lation (see below) we make use of the Healpix Pixelisation 
l|G6rski fc et al.lll999T ). a hierarchical equi-surface and iso- 
latitude pixelisation of the sphere, which was developped to 
analyze polarized CMB type data. 

3 A FULL-SKY MAP MAKER 
3.1 The inverse problem 

Our purpose is now to solve for the non-linear inverse prob- 
lem of recovering the /t(n) map corresponding to a noisy 
incomplete measurement of the 2-D field (g\(n), g2(n)) T of 
the ellipticity and orientation on the sphere (in the local 



tangent plane): 

gfc(n) = .^vL + e fe("). forfc = lor2, (17) 
1 — n(n) 

where n is the sky direction, 7 and k are respectively the 
shear and the convergence, while e is a tensor field of the 
errors which accounts for the measurement noise (includ- 
ing the shot noise induced by the finite number of galaxies 
within that pixel) and model approximations. 



3.1.1 Spherical formulation 

On the sphere, the scalar field n and the tensor field 7 are 
linear functions of the unknown complex field a = Y ■ k 
whose coefficient are the spherical harmonic coefficients of 
k. After discretization and using matrix notation, k and 7 
write 

k = K • a and 7 = G • a , (18) 

where K = Y and G = P Y ■ J, denoting Y the scalar spher- 
ical harmonics and P Y = (_eY,s Y) the parity eigenstates 
based on spin 2 spherical harmonics. These eigenstates are 
defined in such a way that 

71 ± 172 = — ^(a£,m,B ± iai, m ,B)±2 Y e m , 
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Object 


VALUE 


Gain (c-/ADU) 


30.11 


Full well capacity in e- 


300000 


Saturation level (ADU) 


60000 


Read-out noise (e-) 


1.3 


Magnitude zero-point (ADU per second) 


21.254 


Pixel size in arcsec. 


0.2 


Number of microscanning steps 


1 


SB (mag/arcsec 2 ) at 1' from a 0-mag star 


16.0 


Diameter of the primary mirror (in meters) 


8.0 


Obstruction diam. from 2nd mirror in m. 


2.385 


Number of spider arms (0 = none) 


4 


Thickness of the spider arms (in mm) 


5.0 


Pos. angle of the spider pattern 


45.0 


Average wavelength analyzed (microns) 


0.80 


Back, surface brightness (mag/arcsec 2 ) 


21.5 


Nb of stars / D brighter than MAG-LIMITS 


le5 


Slope of differential star counts (dcxp/mag) 


0.3 


Stellar magnitude range allowed 


12.0,19.0 



Table 2. SkyMaker parameters used to generate the VIMOS/ 
VLT images 

so that we have 

= >. , -T.-KT- ae, m ,E+> w + a,e,m,B 

with Wj m = (2Yf_ m ± _2Y^ m )/2. Here J operates on a as 

(J-«W = ( 19 ) 

(J-»)/,m,fl = °- (2°) 

Appendix [B] gives more explicit formulations of the opera- 
tors K and G, using index notation on the sphere. 

3.1.2 Flat sky formulation 

The flat sky limits (corresponding to large ffs) of equations. 
{T8 ]) -(fT9 ]) are (see Appendix EJ: 

J w (1, 0) , and exp(^ • h) , (21) 

while the parity eigenstates read locally, in the fixed copolar 
basis e x ,e y : 

I 2 - I 2 

W + ^-cos(2<^)exp(i£ • n) = - 11 exp(il- h) , 
W"«- i sin(2 4> t ) exp(i £-h) = -i 1 1 exp(i £ ■ n) . (22) 

i ly 

In this limit, the unknowns, a, represent the Fourier coeffi- 
cients of the convergence field, «. Note that our definition of 




Figure 5. Top panel: full-sky view of the mask; Bottom panel: 
a zoom at coordinate (I, b) = (30°, 30°) showing the distribution 
of stellar cuts. This cut corresponds to the inner central region of 
the reconstruction shown in Figure IT3l 

7 and k warrants that they are consistent with the lens equa- 
tion on the tangent plane — solving for k in equation (|18|l 
and plugging the solution into equations (|22|l — which reads 
locally in real space: 

V 2 K (n) = (fl£-flg)7i(fi)+ 2 fl^TiiW. ( 23 ) 

where 71 (x, y) and 72(^,3/) are the two components of the 
E and B modes of the shear field. Also note that thanks to 
equation (|20p the recovered map will not have B modes by 
construction. It can nevertheless be checked that the am- 
plitude of the B modes in the residuals is small compared 
to the amplitude of the signal in the E modes, see Section 
14X31 

3.1.3 Cost function 

The considered problem can be stated as recovering a given 
the data g according to the model in equation (|17|l . In the 
same way as what has been done for deblurring the im- 
ages (section [2} , finding the solution of this inverse prob - 
lem in the Maximum a Po steriori (MAP) (|ThiebautJ (|2005l) ; 
IPichon fc Thiebautl (|l998t )) sense involves minimizing a two- 
term cost function: 

Q(a) = C(a) + ixH{a) , (24) 

with respect to the parameters a. In the right hand side 
of equation (|24|) . the term C(a) enforces agreement of the 
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Figure 6. a zoom of the full-sky recovered k maps of a simulation 2048 Spg 1 * with £ cu t = 722, (top left panel) and 1569, (top right 
panel) (resp. 24, and 78 ngal/D arcmin) at coordinates (<f>,0) = (0,0) (the color table corresponds to a histogram equalization); bottom 
left panel: the corresponding data (the hue color table codes the shear orientation); bottom right panel: the corresponding underlying k 
map. 



model with the data, whereas lZ(a) is a regularization term 
used to enforce our prior knowledge about the sought fields, 
and fi ^ is a Lagrange multiplier used to tune the relative 
importance of the prior with respect to the data. 

For errors with a centered Gaussian distribution, the 
likelihood term writes: 

with e k {hj) = gk(nj) — jk(nj)/[l — n(n)} and W = C _1 
with Cj lihl ,j 2 ,k 2 = {e. kl {n J1 ) e k2 (n j2 )). If the errors are fur- 
ther uncorrelated, the likelihood simplifies to: 



C(a) 



9k(nj 



1 - k(Uj 



(25) 



where the sum is carried over the index j of the sampled 
sky directions fij (so called sky pixels) and index k of the 
two components of, say, the Q and U polarization fields re- 
spectively (see Appendix [B] for an explicit formulation with 
all the relevant indices) and the weights are related to the 
variance of the noise: 



This allow us to account for non uniform noise on the sky 
and also cuts (the galaxy, bright stars, etc.) for which the 
variance can be considered as infinite and thus the corre- 
sponding weights set to zero. Note that setting the weights 
in this statistically consistent way yields no such biases as 
those which would result from interpol ation or inpaintin g 
method s used to replace m issing data (|Pires et al.l 120081 ) , 
see also lAbrial et al.l (|2008h for such implementation in the 
context of CMB experiments). For this recovery problem, 
our prior is that the field k must be as smooth as possible 
in the limit that the model remains compatible with observ- 
ables within the error bars, that is equation (|17[) must be 
valid. To that end, the regularization is written as a penalty 
based on the second order spatial derivatives (Laplacian) 
V 2 k of the field k: 



11(a) = ||V 2 k\\ 



(27) 



w Jyk = Var(e fe (nj))" 



(26) 



Equation (|B12|) in Appendix [Bl gives the expression of V 2 k 
as a function of the unknown a. In order to enforces smooth- 
ness while preserving some sharp features in the k map, 
quadratic and non quadratic norms of the Laplacian have 
been considered for the regularization, see Appendix |B1 
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contrast, in the Born approximation using: 



«(n) = -f2 r 



dz r>(z)£>(z,Z s 



a(z) H 



E{z) V{z s ) 

(28) 

which is valid for sources at a single redshift z s = 1, and 
T>(z) = Hox(z)/c is the adimensional comoving radial co- 
ordinate, hence dT) — dz/E(z). The detailed procedure to 
construct such maps from the simulation using equation (|28[) 
is described in Appe ndix ID II (chosi n g the sampling strat- 
egy) and ID2I and in iTevssier et alj (2008). In practice, a 
set of degraded maps of k was generated from the full res- 
olution, n s idc = 4096 down to n s jd c = 128 in powers of 2, 
together with the corresponding masks (see FigurefS]). Differ- 
ent levels of noise (corresponding to 700 ^ £ C ut < 2500) and 
maps with/without Galactic masks are considered. The cor- 
responding simulations are labeled as 

sian maps are also used, labeled as n . 



SNR 



J FS/GC 



Carte- 



rs 

l°NL/lin 



correspond- 
ing to Cartesian sections of the full-sky maps, where for 
commodity, the experiments involving high resolution where 
calibrated. Here the flag NL/lin refers to whether or not the 
non-linear model is accounted for. 



Figure 7. a zoom on the power spectra of the three reconstruc- 
tions of 2048S'p c s ut for 4 ut = 722, 1083, and 1569 (resp. 24, 44 
and 78 ngal/D arcmin), together with the power spectra of the 
noise. Note that the level of smoothing decreases with increasing 
signal to noise, in parallel to the bias in the corresponding power 
spectrum. 



3.2 Generating the virtual data set 

Let us first describe in turn the simulation used to model 
the full sky k map, and the generation of the corresponding 
map. 



3.2.1 The simulation 

The Horizon 4n (|Tevssier et all (|200St ). iPrunet et ah! 
(|2008D l simulation was used, a ACDM dark matter simu- 
lation using the WMAP 3 cosmogony with a box size of 
2/i~ 1 Gpc on a grid of 4096 3 cells. The 70 billion particles 
were evolve d using the Par ticle Mesh scheme of the RAM- 
SES code (|Tevssierl (|2002l )) on an adaptively refined grid 
(AMR) with about 140 billions cells, reaching a formal res- 
olution of 262144 cells in each direction (roughly 7 kpc/h 
comoving). The simulation covers a sufficiently large vol- 
ume to compute a full-sky convergence map, while resolving 
Milky- Way size halos with more than 100 particles, and ex- 
ploring small scales deeply into the non-linear regime. The 
dark matter distribution in the simulation was integrated in 
a light cone out to redshift 1, around an observer located at 
the center of the simulation box. 



3.2.2 Mock data 

This light cone was then used to calculate the corresponding 
full sky lensing convergence field, which is mapped using the 
Healpix pixelisation scheme with a pixel resolution of A9 ~ 
0.74 arcmin 2 (n s id c = 4096). Specifically, the convergence 
K(n) at the sky coordinate n is computed from the density 



3.2.3 Penalty weight 

In this paper, the weight of the penalty, [i, in equation (|24|l 
is chosen so that the £2 cutoff corresponds to the scale, £ CI n 
at the intersection of the signal and the noise power spectra, 
see e.g. Figure [7] Specifically 

pi oc l/i 2 crit . 

In a more realistic situation, when the power spectrum of the 
signal is unknown, generalized cross validation could be used 
to find this scale. When ti — £2 penalty is implemented (see 
Section |2.1.4| l. the £\ parameter e entering equation (|16[) is 
chosen so that it cuts off the tail of the PDF of the Laplacian 
of the recovered field at the 3-cr level. 

3.3 Optimization &: Performance 

Let us now turn to the optimization procedure and the per- 
formance of the algorithm. 

3. 3. 1 Optimization 

Recall that the procedure assumes here a sampling strat- 
egy, since the noisy g field is given on a pixelisation of the 
sphere. To solve the optimization problem, we used the al- 
gorithm vmlm from OptimPack (|ThiebautI I2002T ) which 
only involves computing the objective function Q(a) and its 
partial derivative with respect to the parameters a. Vmlm 
is an unconstrained version of vmlmb which has been used 
for the deblurring problem and which is described in some 
details in section 12.11 The optimization of equation (|24p is 
carried by computing in turn e quation (11811 and Equation s 
([28)1 and ((B7| using Healpix (|G6rski & et al.I (j 1999T )) 
in OpenMP or MPI. 

3.3.2 Overall Performance 

Each back and forth transform takes respectively 0.1, 0.5, 2, 
8, 32 and 128 seconds on an octo OPTERON for n s id c equal 
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Figure 8. Top panel: a map of the 100 times difference between 
the recovered map with the non-linear model and £\ — £2 penalty, 
and the recovered map without accounting for the non-linearity. 
As expected the difference is largest at high frequencies near the 
cluster and along the filaments. Bottom panel: the power spec- 
trum of the relative difference as a function of £. 




Figure 9. top left panel: a zoom of the original map at coordi- 
nates (l,b) = (0°,0°); top right panel: reconstruction with £\ — £2 
penalty using the non-linear model, bottom left panel: input map 
smoothed at a FHWM of 1.5 pixels, bottom right panel: recon- 
struction with £2 penalty using the non-linear model. The color 
table is linear. The edge-preserving penalty appears qualitatively 
to preserve much better the amplitude and the number of high 
peaks in the re map, as shown quantitatively in Figure ITOl 



cal and topological tracers. We chose a selection of statistical 
tests that are sensitive to different aspects of map-making. 



to 128, 256, 512, 1024, 2048 and 4096, see Table Q] The 
linearized problem without mask converges typically in a 
dozen iterations (which typically only involve a back and 
forth transform, unless the convergence is poor). The lin- 
earized mask problem takes a few hundred iterations, see 
Table |3j and so does the non-linear problem (or the lin- 
earized problem with a non-linear £\ — £2 penalty function) . 



4 VALIDATION AND POST ANALYSIS 

Let us illustrate on a sequence of statistical tests several cru- 
cial features of the ASKI map making algorithm: its ability 
to fill gaps, its ability to preserve the geometry and sharp- 
ness of clusters and maintain the gravitational nature of the 
signal in the presence of masks, and the freedom to choose 
strong/weak prior on the two-points correlation. These prop- 
erties are important in various contexts of the weak lensing 
studies, such as the estimation of cosmological parameters, 
the physics of clusters, the interpretation of tomographic 
data from upcoming surveys, constraining the dark energy 
equation of state through the redshift evolution of statisti- 



4.1 One point statistics 

4.. 1.1 Cluster counts 

One of the main assets of high resolution full-sky lensing 
maps is to probe multiple scales: it then becomes possible 
to sample the non linear transition scale and, e.g. study the 
shape of clusters. Figure [9] illustrates this feature while dis- 
playing the result of the inversion with £2 and £\ — £2 penal- 
ties. For this experiment, a Cartesian subset at galactic co- 
ordinates (l,b) — (0°,0°) was extracted. The corresponding 
non-linear shear field g was generated via Fourier transform, 
and noised with a white additive noise of SNR of 1. This 
set was then inverted while assuming £2 (bottom right) and 
£1 — £2 (top right) penalties. The choice for the two penalty 
weights, /j, and e was made on the basis of least square resid- 
ual in the inverted k maps. The improvement of £1 — £2 over 
£2 penalty is significant. This statement is made more quan- 
titative in Figure \W\ which displays the PDF of the peaks 
within that image for the initial map (top left panel of Fig- 
ure [9]) computed following the peak patch prescription de- 
scribed in Section 14.3.11 The agreement between the input 
and the recovered distribution is significantly enhanced by 
the optimal £\ — £2 (top right) penalty. 
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Figure 11. left panel: skewness, 53 and kurtosis, S4, as a function of scale (using sharp top hat filtering) for the model (plane line) and 
the recovered k maps of a simulation 2048 ^ps" (dotted, dot-dashed, dot- dot- dashed line for £ cu t = 722, 1083, and 1569, resp. 24, 44 and 
78 ngal/D arcmin); right panel: same as top panel, but for the 2048 Sqq set. Note that the kurtosis of the cut is significantly different at 
small I. 



Point Source Distribution with/without L1-L2 penalty 




K 

Figure 10. The PDF of point source as defined in Section 11.3.1 1 
corresponding to the maps of Figure [9] recovered with i\ — £2 
penalty, and £2 penalty respectively. The improvement with an 
edge-preserving penalty is significant. 

4.1.2 Skewness and Kurtosis 

The simplest statistics to explore the non linear transition 
is the skewness, 53 and the kurtosis, S4 of the PDF of 
the recovered maps. Furthermore, it has been shown that 
these parameters provide a powerful tool to measure the un- 



derlying cosmological param eters (|Bernardeau et all (|l997t ) ; 
iTakada fc jaml (20021 , 120041 ')'). Figure ITT1 displays the evolu- 
tion of these numbers as a function of scale in the initial and 
recovered maps, with and without galactic masking. The 
top hat filter used here is of width [2% 2 8+1 ], while the har- 
monic number of each band is the mean of its boundary: 
i = (2 1 + 2 l+1 )/2. The recovery of skewness and kurtosis 
is good in the case of unmasked data. Of course it degrades 
with the scale as we reach £ cu t- Using the reconstructed map 
is not the optimal way of measuring the 3 and 4 point func- 
tions at small scale. However, an optimal dedicated estima- 
tor can be built upon the same regularization technique. The 
masked case is not as good. There, a dedicated estimator, 
acting only on small, clean, pieces of the sky will probably 
yield better results. 

4-1.3 Accounting for a non-linear model 

Figure [8] shows the effect of accounting for the non linear- 
ity in equation (|17|1 . Here a set of Cartesian simulations 
is used 256C' NL / lin . This map represents (a 100 times) the 
difference between the recovered map while accounting for 
1 — k in equation (|17|l in the inversion, and the recovered 
map while neglecting this factor. The difference is small in 
amplitude, but shows as expected the strongest bias near 
the clusters and the filaments, where k is largest. The bot- 
tom panel represents the corresponding relative power spec- 
trum, Ci [NL — linJ/CV [input] as a function of I. Again the 
larger discrepancy occurs at higher I, corresponding to the 
sharp peaks at the positions of the clusters. Hence the 
non-linearity should be accounted for in th e mode l if th e 
shape of the clust e r is an issue ( s ee als o IWhitd (|2005h : 
iDodelson fc Zhand (120051 ); IShapirol (|2009h V For all practi- 
cal purposes, we have therefore demonstrated that at scales 
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below imax < 4096 solving the linearized problem is de facto 
equivalent to the general non-linear problem when k is ne- 
glected at the denominator in equation (|17[1 , 



4.2 Two points statistics 

Since ASKI was constructed to provide the optimal map 
given the measured shear, we do not expect that it will yield 
the optimal estimator for non-linear functions of these maps, 
such as the powerspectrum, bispectrum etc.. Nevertheless it 
is of interest to compare the two point statistics of input 
and recovered maps, to see how Aski deals with masks and 
how it affects the occurrence of spurious B modes. 



4-2.1 Optimal Wiener filtering 

Throughout this paper the prior Ci = £~ 1 (£+ ("Lapla- 
cian prior") is used in equation (|B11[) . Let us briefly inves- 
tigate how a customized (Wiener) prior for Ci changes the 
reconstruction at small scales. Figure [P2l shows that the cor- 
responding power spectra of the reconstructed kappa maps, 
as expected, differ mostly for scales where the signal to noise 
is smaller than one. However, when the smoothing (Lapla- 
cian) prior amplitude (see equation I27|l is tuned to minimize 
the reconstruction error as in the figure, the power spectra of 
the reconstructed maps for the two different priors (Lapla- 
cian and Wiener) are quite similar. 

It is interesting to note that the power of the reconstructed 
map with the Wiener prior (light brown line in Figure I12[) 
is systematically biased low as compared to the input power 
spectrum. This reflects the fact that an optimal (minimum 
variance) estimation of the power spectrum is not equivalent 
to a power spectrum estimation on an optimal (minimum 
variance) reconstructed map. However, in the simple case 
where we have noisy data without masks, the bias of the 
power spectrum of the minimum variance map reconstruc- 
tion is known, it is simply given by Ct/(Ce + Ne) where Ce 
is the power spectrum of the underlying kappa map (with- 
out noise), and Nt is the noise power spectrum in "kappa" 
space, which is given approximately in our case by a 2 Q, P i X , 
where a 2 is the noise variance per pixel in the shear field, 
and flpix is the solid angle of a pixel. 

In Figure [T^l the noise level is shown by the horizontal dark 
blue line. One can see in particular in the figure, that when 
the model power spectrum (without noise) crosses the noise 
power spectrum, the power spectrum of the minimum vari- 
ance map (golden line) is lower than that of the model by 
a factor of 2, as expected from the considerations above, 
even in the presence of masks. Thus an approximate, but 
simple way to get an unbiased estimate of the kappa power 
spectrum is to correct the minimum variance map power 
spectrum by the ratio Ce/(Ct + Ni). Note however that a 
true minimum variance power spectrum estimation of the 
kappa field is not the aim of the present method (see e.g. 
[Pent (|2003T > for the flat-sky case). 

Nevertheless elsewhere in this paper, a smoothing prior 
which is not customized to the specific problem is preferred. 
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Figure 12. The power spectrum of the input and recovered re, 
(with smoothing and Ce prior, see equation l)B8|l ) as a function 
of I, together with the power spectrum of the noise and the noisy 
equivalent re using a simulation, si^Spg 4 (i.e ngal/Darcmin = 5). 
Note that the recovered power spectrum departs from the power 
spectrum of the input field roughly at the cutoff frequency when 
a quadratic smoothing penalty is applied. 
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Table 3. same as Table [Tj with Galactic masks. 



4-2.2 Filling gaps within masks 

Let us first compare visually the recovered map to the input 
map. Figure [13] illustrates a feature of the penalized recon- 
struction: it interpolates quite well and provides means to 
fill the gaps corresponding to the galactic cuts. For a more 
quantitative comparison we also plot on this figure the ridges 
of both fields (using the skeleton, see below), which match 
very well up to the very edge of the mask. The smoothing 
penalty also induces a level of extrapolation, best seen in 
the residuals, see Figure 1181 The masking (or more gener- 
ally, non uniform weights, ua) nevertheless biases the re- 
constructed map, as seen on Figure [TT1 and [T4l Note finally 
that when masks are accounted for, it is straightforward to 
correct for them when computing the powerspectrum as the 
harmonic transform of the autocorrelation, which in turn is 
deri ved by correcti n g for t he autocorrelation of the masks: 
(see ISzapudi et all (|200lD : iHivon et al l (|2002r i; IChon et all 
|2004[ ) for details). When seeking the three-point correla- 
tions one could also proceed accordingly, and divide by the 
three-point correlation of the mask. Indeed a three-point 
reduced correlation is simply one plus the excess probabil- 
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Input k: l=30,b=30 Recovered map: l cu t=221 2, 1=30, b=30 Skeleton of k with masks 
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Figure 13. left panel: the initial k map in the region with bright stars masking shown in Figure [5] at coordinates (I, b) = (30°, 30°); 
middle panel: the corresponding recovered k map of a simulation 2048 Sqc 2 Note that the gaps have been nicely filled up to the very 
edge of the mask; right panel: the corresponding two skeletons (color coded by k in purple: input skeleton; in orange: recovered skeleton) 
for the inner region (marked as a square on the middle panel), when masking is present. Note the clear gradient away from the mask in 
the quality of the match between the two skeletons; Recall that most of this field is partially shielded by stars, as seen in Figure [5] 



Powerspecirum: linear model, with mask 
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Figure 14. the power spectra of three high resolution reconstruc- 
tions corresponding to Figure fl3l for £ cu t = 796,1368, and 2212 
(resp. 28, 63 and 130 ngal/D arcmin corresponding to a low, in- 
termediate and high end values) together with the power spectra 
of the noise. Note that the the recovered power spectrum has ex- 
tra power at large scales and less power at intermediate scales, 
an artifact of the mask which can be corrected for by accounting 
for the prior knowledge of the auto-correlation of the mask. 



ity of finding triplets, which in turn is computed by count- 
ing the number of found triplets and dividing by the ex- 
pected number of such t riplets given the shape of the mask 
l|Chen fc Szapudil (|2005h . This also applies if the mask is 
grey 



4-2.3 Residual B modes 

Let us investigate the effect of leaking of B modes with the 
following experiment: the noise in the transform of the B 
channel is boosted by some fixed amount over a map which 
has Galactic cuts. This corresponds to the case where the B 
is significantly larger than the noise, yet uncorrelated with 
the E mode, corresponding to e.g. a systematic bias in the 
ellipticity extraction for example. It is expected that, due 
to masks, this B mode will leak in E. An example of such 
leak, for B modes as large or up to 32 times larger than 
the noise, is shown in Figure [18] The power spectrum of 
the residuals in the corresponding k map is computed while 
masking in the residual the exact regions corresponding to 
the cuts. When this boost is zero, (bottom curve in Figure 
I17p the power spectrum of these residuals is flat and corre- 
sponds to the noise powerspectrum. In contrast, the stronger 
the boost the larger the scale below which this power spec- 
trum is colored. Note that it was checked that, as expected, 
these coherent residuals disappear completely if the galactic 
cuts are ignored. It would also be interesting to compare the 
distribution of the shape of dark matter in input/recovered 
clusters. 

Finally, note that Appendix IC2I discusses briefly the ef- 
fect of noise in powerspectrum estimation. 

4.3 Alternative statistics: critical sets 

Let us close this section with a quantitative comparison of 
the input and the recovered map using more exotic probes to 
estimate the quality of the reconstruction, and the prospect 
it offers for dark energy measurements. Indeed, the predic- 
tions of the perturbative hierarchical clustering model are 
often given through the hierarchy of the differences between 
the moments to their Gaussian limit. Yet higher order mo- 
ments are generally difficult to test directly in real-life ob- 
servation s, due to their sensitiv ity to very rare events. As 
argued in lPogosvan et all (|2009l ) the geometrical analysis of 
the critical sets in the field (extrema counts, Genus, criti- 
cal lines etc..) may provide more robust measures of non- 
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Gaussianity, and is becoming elsewhere an act ive field of 
investigation l|Gott et al.ll2009l ; IPark et al.ll2005l ). 

4-3.1 Peak patch counts and area 

Even though many tools are available to identify peaks 
within the reconstructed map, let us validate here our recon- 
struction using a segmentation of both the initial and the 
recovered maps using peak patches on the sphere, which are 
a segmentation of the map based on the attr action patches 
of the k map when following its gradient (see lSousbie et all 
(|2008h V Within each peak patch (see Figure iBlj) . the bright- 
est pixel is assigned a mass corresponding to the enclosed 
mass within the peak patch. This quantity is gravitation- 
ally motivated (as the patch corresponds to the attraction 
region of the cluster) and is both robust (as the geometry of 
the patch only depends on the imposed smoothing length, 
which in turn is fixed by the resolution of the survey) and 
sensitive to small features in the map; it is therefore a good 
indicator of the quality of the reconstruction. Figure [Tol (left 
panel) displays the corresponding PDFs before and after re- 
construction. As expected, the recovered point source PDF 
has a shifted mode and is less skewed than the original dis- 
tribution. This trend decreases with increasing SNR. For 
realistic galaxy counts of 40 ngal/D arcmin, the agreement 
between the input and the recovered PDF is fairly good, 
and the corresponding residual bias can be modeled (as the 
reconstruction is essentially a smoothing of the underlying 
map). This could lead to interesting constraints on f2 m and 
D(z) when used in conjunction with weak lensing tomog- 
raphy in order to probe its redshift evolution. The right 
panel of Figure [15] focuses on a different quantity, the area 
of the patches, which when compared to the area of the cor- 
responding void patches, could also be used as a measure of 
the gravitationally induced non gaussianities, together with 
their shape (higher moments of k within a patch). Again, the 
reconstruction seems to recover this distribution well enough 
to suggest that such a tool could be used in the future to 
study the cosmic evolution of the projected web. 

4-3.2 Topology & geometry: critical lines 

Let us now compare the shape of the recovered map to the 
initial map from the point of view of its critical lines. For this 
purpose, let us use h e re the skeleton as a geometric probe 
|Novikov et al.l (120061 ) ; ISousbie et al.l (|2007h ). It is defined in 
2D as the boundary of the void patches, which in turn are a 
segmentation of the map based on the valleys of the k map 
(corresponding to the peak patches (defined above) of minus 
the field) . The skeleton of the initial field and the recovered 
fields for simulation 2048S , qc : is computed, and represented 
in Figure [13] The recovered skeletons are qualitatively fairly 
close to the original skeleton, which demonstrates that the 
local topology and geometry of the field is well recovered. 
Let us make this comparison more quantitative. The dif- 
ferential length per unit area of the recovered field (the set 
204s5'p o s ut with 4ut = 722, 1083, and 1569 as labeled^ over 
the initial k map (thin line) as a function of density thresh- 
old is also shown in Figure \W\ while Figure [6] shows the 

9 Note that ngal/Darcmin = 40(Z C ut /1000) 1 ' 5 . 
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Figure 16. The input skeleton differential length (a tracer of 
f! m ) with its recovered counterparts as a function of the nor- 
malized k contrast, v = (k — R)/a K for the set 2048 •S'pg 1 ' with 
£ cut = 722, 1083, and 1569 (resp. 24, 44 and 78 ngal/D arcmin) . 
Here the PDF of the normalized n contrast was subtracted to the 
differential length for clarity. As expected, the agreement is best 
at large convergence. This figure is complementary to Figure [13] 
which shows that the geometry of the field is well preserved on 
average. 



corresponding maps for similar runs, together with a map 
of the orientation of the g field. The agreement increases at 
larger density thresholds, which s ugg ests that the topology 
of dens e regions is w e ll reco vered 10 !. The total length was 
shown (|Sousbie et al.l (|2008l )) to trace well the underlying 
shape parameter of the powerspectrum and has been used 
in 3D to constaint th e dark matter content of the universe 
(|Sousbie et all (|2006h ). As shown in iPogosvan et al.l (|2009l ) 
this would work for two dimensional maps and could there- 
fore be used with k maps such as those reconstructed via the 
present method. The redshift evolution of this differential 
count, when tomographic data is available, could comple- 
ment e.g. Genus measurements as means of constraining the 
dark energy equation of state in a manner which could be 
more robust than direct cumulant estimation. Eventually, 
the skeleton could also be used to characterize the connec- 
tivity of clusters (i.e. the number of connected projected 
filaments), as it will als o depend on t h e cosm ic dark energy 
content of the universe l|Pichon fc al.l (|2009l )). 

This rapid review has shown that, depending on the 
final objective (cosmological parameters, cross correlation 
with other maps etc.), a variety of estimators can be ex- 
tracted from the recovered maps. Aski was shown to per- 
form rather well with respect to these estimators. Defining 
the best combination of these estimators, - and the opti- 



In fact the relative distance between the recovered and the in- 
put skelet on could also be used as an alternative to the differential 
length see ICaucci et alj d2008h . 
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Figure 15. left panel: the PDFs of re max at point sources before and after reconstruction of set of simulations 2048 ^pg" (dashed, dotted, 
dot-dashed line for £ cu t = 722, 1083, and 1569 (resp. 24, 44 and 78 ngal/TJ arcmin); right panel: the PDFs of the area of peak patches 
(see Figure lB"lT l before and after reconstruction for the same set of simulations. Note that, as expected, the recovered distribution of 
peaks is less skewed than the original, whereas conversely, the PDF of the area of the peak patches for the low SNR reconstruction is 
more skewed towards larger patches. 



mal penalty associated - will be one of the key topic lensing 
research for the coming years. 



5 CONCLUSION & DISCUSSION 

This paper sketched possible solutions to issues that a full- 
sky weak-lensing pipeline will have to address, and presented 
an inverse method implementing the debluring of the image 
and the map making step. 

Weak lensing surveys require measuring statistical dis- 
tributions of the morphological parameters (ellipticity, ori- 
entation, ...) of a very large number of galaxies. This pa- 
per demonstrated that these parameters can be measured 
with a better accuracy and strongly reduced bias if the deep 
sky images are properly deblurred prior to the shape mea- 
surements. Using a relative figure of merit (the recovered 
SExtractor ellipticity) we have shown that this deblur- 
ring could in crowded fields improve more than tenfold the 
accuracy of the recovered ellipticities. This deblurring is crit- 
ical in crowded regions, where the overlapping of stars and 
galaxies otherwise prevents accurate morphological estima- 
tion. Henceforth dealing with such regions is important for 
a full-sky survey. Since such surveys will require the pro- 
cessing of a great number of large images, the calibration 
of these techniques is automated on the images themselves 
via cross validation after identification and removal of the 
stars within the field (see Figure 1). In particular the level 
of regularization, fi and the l\ — £2 threshold are automat- 
ically tuned in order to deal with the noise level and the 
dynamics of the raw images. The gap-filling interpolation 
feature of the inversion would apply even more efficiently 
in this regime than in the map reconstruction regime de- 
scribed in Section [3] The algorithm described here scales 



Residual with boosted B mode 




Figure 17. Power spectrum of the mask weighted residual er- 
ror on k as a function of the harmonic number, i. The different 
curves correspond to boost of the B modes of increasing relative 
strength. The low order modes are polluted by leaks from the 
masks (see also Figure \18\ : here ( cu t = 752 (60 ngal/D arcmin). 

well since it only relies on DFTs: hence it could be applied 
to very large images such as those produced by modern sur- 
veys. Aski uses the efficient variable metric limited mem- 
ory algorithm OptimPack, which allows both optimizations 
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to scale to high resolutions. The debhirring is implemented 
on Cartesian maps as large as 16 384 2 pixels. Generalized 
cross validation was shown to yield a quantitative threshold 
in order to remove accurately the point sources within the 
field, hence imposing the optimal level of smoothing for the 
galaxies only. In this paper, the focus was put on blurred 
8-meter ground-based observations, but the implementation 
for EUCLID-like space missions should be straightforward. 
The above described improvements could clearly be repro- 
duced if alternative state of the art shear estimators were to 
be used (as compared by the SHear Testing Program). 

This paper also demonstrated that optimization in the 
context of Maximum A Posteriori provided a consistent 
framework for the optimal reconstruction of k maps on the 
sphere. The main asset of the Aski algorithm is that the 
penalty can be applied in model space, while the optimiza- 
tion iterates back and forth between data space and model 
space. This freedom allows it to deal simultaneously with 
masks (in data space) and edge preserving penalties. Pro- 
viding k maps is critical both in its own right, as it maps 
the dark matter distribution of our universe and gives access 
to the underlying powerspectrum at large scale. Such maps 
are also interesting when cross-correlated with other surveys 
(optical surveys, CMB maps, lensing reconstruction and dis- 
tribution of SZ clusters from the Planck mission, redshift 
evolution of X-ray sources counts etc..) in order to explore 
the evolution of the large-scale structure, and in the case 
of the surveys mapping the baryonic matter, to better un- 
derstand biasing as a function of scale. Finally, though not 
optimally, it can be used to compute second and higher or- 
der statistics, and noticeably the three-point statistics, the 
Genus, cluster counts or the skeleton, which may constrain 
more efficiently the dark energy equation of state, as they 
are less sensitive to rare events. It should be stressed once 
more that while the reconstructed k maps yield biased esti- 
mates of the power-spectrum and higher order statistics, the 
technique described in this paper can be adapted to build 
dedicated optimal estimators for each of those observables. 

Section [4] demonstrated the quality and limitations of 
the reconstruction using various statistical tools on a full- 
sky simulation of g with resolutions of up to 12 x 4096 2 = 
201 326 592 pixels thanks again to the efficiency of Optim- 
PACK. In particular, it identified point sources of the fields, 
analyzed their PDF and showed that £i — £2 penalty was 
critical at small scales. It also investigated the effect of leak- 
age of B modes when Galactic cuts are present. It presented 
a method to probe the topology and geometry of a field on 
the sphere, the peak patches and the skeleton, and applied 
it to compare the recovered field to the initial field. Such 
tools allow us to quantify the differences between the two 
maps and act as an efficient source segmentation algorithm. 
Indeed, the degeneracy between the cosmological parame- 
ters (r2M,o"g) is for instance best lifted with cluster counts. 
They may also turn out to be of importance when probing 
the dark energy equation of state as they are less sensitive 
to rare events. The Cartesian dual formulation of ASKI was 
also implemented and may prove useful for surveys where 
sky coverage is sufficiently small. 

In short, Aski accounts for the possible building blocks 
that a full-scale pipeline aiming at sampling the dark matter 
distribution over the whole sky should provide. Specifically 
it allows for (i) automatically deblurring very large images 




Figure 18. shows an example of full-sky leak of the B modes 
when masks are accounted for; top panel: the residuals corre- 
sponding to ag = a + a: the inner box corresponds to a zoom 
near the edge of the galactic cut at (b, I) = (30, 20); bottom panel: 
same residual and box for ug = a + 32ct. Note that for the latter 
case, the extend of the leakage is much larger and coherent. 

using non-parametric self-calibrated edge-preserving l\ — £2 
deconvolution with positivity; (ii) carrying the large non- 
linear inverse problem of reconstructing the convergence k 
from the shear g using equation (JTTJl : the back and forth 
iterations between model and data are consistent with con- 
straints in both spaces, and allow for an accurate recovery of 
cluster profiles and shapes; (iii) non-uniform weighting and 
masking: consistent with realistic Galactic cuts (and bright 
stars masking) and non-uniform sampling of the different 
regions of the sky, dealing transparently with the issue of 
the boundary; (iv) edge-preserving i\ — £2 penalty yielding 
quasi point-like cluster reconstruction. Finally (v) it intro- 
duced peak patches and the skeleton on the sphere, together 
with its statistics. 

Possible improvements/investigation beyond the scope 
of this paper involve: (i) comparing the absolute gain in 
shear est imation using alter native tools to Sextractor 
(such as iMassev et alj (|2007)) with more realistics galac- 
tic shapes; (ii) deblurring the images with a variable PSF 
within the field; (iii) building optimal estimators for the 
power spectrum Cf, or the asymmetry 5*3 (a possible op- 
tion would be to rely on perturbation theory, and invert 
the non- linear problem for both Cj> and 53); (iv) inverting 
for 7 and k simultaneously and checking a posteriori the 
amplitude of the B modes (an alternative to the model de- 
scribed in equation (|19|l ; the issue of unicity of the solution 
will be a challenge); (v) carrying the deprojection while as- 
suming prior knowledge of a complete distribution of source 
planes in equation (|28|l (the corresponding inverse problem 
remains linear, with an effective kernel which depends on 
the optical configuration and the distribution of galaxies as 
a function of redshift); (vi) moving away from the Born ap- 
proximation, which involves solving Poisson's equation for 
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each slice, and ray-tracing back to the source while solving 
for the lens equation though all the slices; (vii) implementing 
a more realistic noise modeling (which amounts to changing 
the cost function, equation (|24l) ); (viii) studying the shape 
of dark matter distribution in clusters and groups: typically 
this would also involve cross-correlating the corresponding 
distribution with the light at various wavelengths, (ix) defin- 
ing the post analysis which most sensitive to dark energy, 
given the feature of the surveys to come and finally (x) prop- 
agating the analysis up to the cosmic figure of merit for the 
dark energy parameters. 
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APPENDIX A: EFFICIENT STAR REMOVAL 

We have observed that for realistic deep field images, gen- 
eralized cross validation (GCV) yields an hyper-parameter 
value which is relevant to regularize the higher part of the 
dynamic (mainly due to stars, i.e. point-like objects which 
concentrate their luminous energy in a very small area) but 
which is much too low to regularize the lower parts ot the 
dynamic where galaxies remain. Indeed, when dealing with 
images with a large dynamical range, GCV yields a value of 
the regularization level n which is necessarily a compromise 
between not smoothing too much the sharp features and suf- 
ficient smoothing of low contrasted structures to avoid noise 
amplification. The solution to the problem of underestimat- 
ing the regularization weight can be solved by applying the 
GCV method onto the image with no stars. We want to find 
structures of known shape s(x) but unknown position and 
intensity in the image y. In our case, s(x) is the PSF since 
we want to detect stars. This reasoning could however be 
generalized to other kind of objects. If a single object of 
this shape is present in the image, this could be achieved by 
considering the objective function: 



</>fuii(a, t) — Wk [ a s ( x k — t) 



Vk\ 



structures (or with other fainter structures), a better strat- 
egy is to limit the local fit to a small region of interest (ROI) 
around the structure. This is achieved by minimizing: 

<f>(a, t) = ^w k r(x k - t) [a s(x k - t) - y k f , 

k 

where r(Sx) is equal to 1 within the region of interest (ROI) 
and equal to outside the ROI. Minimization of 4>{a, t) w.r.t. 
a yields the best intensity for a local fit around t: 

d(f> _ ^, _^ k w k r(x k — t) s(xk — t)y k 

da Y^k Wk r ( Xk ~ *) s ( Xk ~ t) 2 

Inserting a* in the objective function yields: 

<f(t) 4 <f>{a,t)\ a = a * > 

= 2_^w k r(x k -t)y k 



J2k Wk r ( Xk ~ *) S ( Xk - t) 2 

Since r(5x) 2 = r(Sx), by defining snoi(Sx) = r(8x) s(8x), 
the local criterion and local best intensity can be rewritten 



0*(t) 



«*(*) = 



^r(x k -t)w k y 2 k 



(J2k s noi(xk - t) w k y k ) 



X] fc SROi(a;fe - t) 2 w k 
J2k s noi(x k - t) w k y k 



to be minimized w.r.t. the weight a and the offset t, here 
a 2-D vector. In fact, since y may be crowded with similar 



J2k s ROi{x k ~ t) 2 w k 

These parameters can be computed for all shifts by an in- 
teger number of pixels by means of FFT's (cross-correlation 
product). Unfortunately, the overall minimum of <j>*(£) is not 
the best choice for removing the brightest structures since 
there is no warranty that this minimum corresponds to a 
bright object. It is better to select the the location which 
yields the brightest structure, i.e. the maximum of a*(t). 
After removal of the contribution a*(t*) s(x — t*) from the 
data, this technique can be repeated to detect the second 
brightest source, and so on. The co rresponding algorithm is 
very similar to the CLEAN method (|Hogbomlll9"74T ISchwarzl 
Il978f) with the further refinement of accounting for non- 
stationary noise and missing d ata. It has been sh own that 
it achieves sub-pixel precision (|Soulez et al.ll2007l ) and that 
it could be used to detect (and remove) out of field sources 
l|Soulez et al.ll2007t ). 



APPENDIX B: MODEL ON THE SPHERE 

Let us describe in more details the model used for the in- 
version of Section [3. II 



Bl Discretization and Sampling 

After discretization and using explicit indices, the model in 
equation (|17[) writes: 

7j,fe , 
Sj.fe = -r^ — + e j,k , 

1 Kj 

where the index j runs over the sky coordinates Aj = 
(xj,yj), index k corresponds to the two components U and 
Q of the polarization, whereas I and m are the harmonic 
indices and p refers to the two components of the spinned 
2-harmonic. In words, the discretization yields: 

9j,k = 9k(nj), 7j, & =7 fc (nj), K=K(nj), e j>k = e h {nj) . 
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Figure Bl. Peak patch of the recovered n map. The inner box zooms the central region. The color coding corresponds loosely to the 
density of the different peak patches. The PDF of the area of these patches is described in Figure IT5l while the maxima mentioned in 
this figure are found within each patch. 



Here the fields k and 7 are linear functions of the complex 
field a of the spherical harmonic coefficients of n. Using the 
matrix notation of the paper, k and 7 write: 

k = K • a , 7 = G ■ a , 

where K = Y and G = P Y- J; with explicit index notations: 

and 

7j,fc = Gj,k,i,m ae,m = ^2 jYj,k,e,™, P (J • a) e m p , 

JZ,rn £,m,p 

To get the detailed expression of the operator J we start 
from the relationship between the lensing potential, the con- 
vergence and the shear fields on the sphere. To do this, we 
need first to define the null diad, based on the polar coordi- 
nates unit vectors: 



m 



(e e + i&j 
V2 



(Bl) 



Given this diad, the lensing potential, convergence and shear 
are related through: 

ViVjtl? = ngij + (71 + i72)(m+ ® m+)ij 
+ (71 — J72)(m_ (g> m_)y , 

where <?y is the spherical metric tensor, and V the spheri- 
cal covariant derivative. Now, using the following expression 
of the second covariant derivative of a scalar spherical har- 
monic: 



ViV,-Y. 



1 (1 + 2)! 



2 V 0-2)! 



[ 2Y £jm (m + ® m+) 



+ - 2 Y tm (m-®m-) Ye t7n gij 

J ij I 

we can relate the convergence and shear fields to the spher- 



ical harmonic coefficients 4v m of the lensing potential: 

= -Ek l+1 )*'.™Y^(ft), (B2) 



(7i±*72)(n) = J2o 



1 1 (1 + 2)} 

s 2Va-2)i 



$e, m ±2Ye irn (h), (B3) 



(n), (B4) 



where the last equality defines the E and B modes coeffi- 
cients. Relating the latter coefficients to the spherical har- 
monics decomposition of k written above, we get the follow- 
ing expression for the J operator coefficients: 



(J • a) e 
(J • a) t 



(* + 2)(l-l) 
(£ + 1)1 



= 0, 



(B5) 
(B6) 



where ag^ m are the spherical harmonics coefficients of the 
convergence field. 



B2 Likelihood 

The data related term in the cost function is 

2 



1 — AC, 



9j,k 



The gradient of this term is needed to find the solution of 
the inverse problem: 



dC(a) 



where 



3,k J j (I-"!) 



1j,k 

r 3 , k = Wj, h [ — - gj, k 
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are the weighted residuals, and where 



H 



£,m,j,k — 



(l + 2)(l- l) 
{£ + !)£ 



' P X 1,771,1, j,k ■ 



(B7) 



B3 Regularization 

The aim of the regularization is to avoid ill-conditioning and 
noise amplification in the inversion. Following a Bayesian 
prescription, this can be achieved by requiring the field k 
to obey some known a priori statistics, or while assuming a 
roughness penalty for 1Z. 



B3. 1 Wiener filter and £2 penalty 

Assuming the field k has Gaussian distribution with mean 
k — (k) and covariance C K = ((k — «)•(« — k) t ), the prior 
penalty should write: 

-1 



/'• 



1 and 1Z ■ 



K) 



K) 



For a field with zero mean (k. = 0) and stationary isotropic 
statistics, the regularization can be expressed in terms of the 
harmonic coefficients: 

2 



71(a) 



-1/2. J| 2 _ Y- E, 



2^ 



c, 



with 



Ce = (la* , 



(B8) 



(B9) 



where the angular brackets denote here the expected value 
taken over the index m of the harmonic coefficients. The 
gradient of the stationnary isotropic Gaussian regularization 
in equation (IB8|l is: 

dTZ(a) _ ae, m 
dae.m Ce 

Note that the regularization in equation (|B8|) with a known 
power spectrum Ce for the field n yields the so-called Wiener 
filter. 

When the power spectrum of k is not exactly known, a 
quadratic prior can alternatively be used. For instance: 

2 



n{a) 



R" 



-1/2 



E 



R, 



(BIO) 



In our framework, effective regularization is achieved by re- 
quiring the field n to be somewhat smooth. In practice, this 
is obtained by requiring Re to be a positive non-decreasing 
function of the index I. Note that, from a Bayesian view- 
point, the regularization in equation (|BI Oj) corresponds to 
the prior that k is a stationary isotropic centered Gaussian 
field with mean power spectrum Ce = Re, which is similar 
to the Wiener filter except that the exact statistics is not 
known in advance (because some parameters of the regular- 
ization have to be tuned; for instance, /x need not be equal 
to one). The gradient of 1Z in equation (|BI0[) reads: 



dTZ(a) 



= 2 



ae,n 



The quadratic prior in equation (|B10|l can be expressed in 
terms of k: 

Tl = ||R~ 1/2 ■ all 2 = IID • k|| 2 , 



where D = R~ 1//2 • Y* is some finite difference operator 
which gives an estimate of the local fluctuation of the field, 
and Y* is the (pseudo-) inverse of the scalar spherical har- 
monics matrix. In our framework, we choose to measure the 
amplitude of the local fluctuations of the field n by its Lapla- 
cian V 2 k and to express the regularization penalty as: 

(Bll) 



ft(a)=£>((V 



where the cost function <j>{r) is an increasing function of \r\. 
When 4>{r) = r 2 , our regulariztion is a quadratic penalty 
similar to equation (|BI0[) . Using matrix notation, the Lapla- 
cian of the field n write : 



V k = YL 



-1/2 



with 



-1/2 



17' 



where Le = £~ 2 (£+ l) -2 . In order to perform the minimiza- 
tion, the gradient of the regularization must be computed. 
By the chain rule: 



dTZ(a) 

dae.m 



dae-. 



= E 



V 2 «) 



(B12) 



where <j>'{r) is the derivative of 4>{r). 



B3.2 £2 - £1 penalty 

As for the image restoration, quadratic regularization yields 
spuriours ripples in the regularized n map. To avoid them, 
we propose to use a £2 — £1 cost function <j> applied to the 
Laplacian of n. The details of the £2 — £1 cost function are 



discussed in section l2.1.4I Taking IZ(a) = • 
with cj> given in equation (IB) , yields: 



dTZ(a) 

dae.m 



= E 



2e (V 2 k) . 9(V 2 /t) 



£+ (V 2 /t) 



2e E- 



■ £,m,j 



(!H 

£+ |(V 2 /t) 



In practice, we use GCV to set the level of the regularization, 
possibly after cluster removal (as explained in Appendix |A"| 
and the £\ — £2 threshold is set to be e — a a where a ~ 2 — 3 
and a is the standard deviation of the histogram of spatial 
finite differences. 



APPENDIX C: FROM THE SPHERE TO THE 
PLANE 



In section lB.f ,2l we sketched the correspondence between the 
fullsky and the flat sky approximation of the lens equation. 
Let us derive it here precisely and use it to investigate the 
effect of shot noise in the estimation of k. 



CI Derivation 

Following closely iHul (|2000l ). let us start with a scalar field 
on the sphere, and its decomposition on the usual spherical 
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harmonics: 



x(h) = Y, Xe '™ Y 



and let us define 



-Ity 



2f- 



t,mt 



i4> l 



(CI) 



(C2) 



together with the inverse relation 
Xi m — 



21+1 



4tt 



2n 



where <f>e is the polar angle of the 1 vector in Fourier space. 
Let us show that X(l) corresponds to the Fourier decompo- 
sition of the field in the flat-sky limit (small angles near the 
pole). Indeed, taking the asymptotic behavior of the spher- 
ical harmonics 



together with the plane-wave expansion in terms of Bessel 
functions 



e ii.ft = Y^i m J m (£9)e xrn{4 - M f 

m 

We get from equation (|C1[) 

B J 



x(iy ln . 



irrupi 



(277) 

For a spin-2 field, let us proceed likewise. We start from the 
all-sky definition of a spin-2 tensor field, and its decomposi- 
tion in spin-2 spherical harmonics: 



±X(h) = 2_. ±Xt. m ±iY t. m , 



(C3) 



where ±X(n) is defined in the spherical tangent coordinates 
eg,e^. We define, as in equation ()C2[) , the Fourier modes of 
the components of the spin-2 field as ±X(l). We have in the 
flat-sky limit the following asymptotic form for the spin-2 
spherical harmonics: 



± 2 Y l , m « ^2 e ± id v) Y *,™ . (C4) 

Plugging equation ()C3|I into equation ()C4[) yields: 
±A(n) 



2—i 9-jr / O.TT ^ ' 



ri 2 l 

(277)= 



±X(l)c 



±2i(0 f - 



Redefining the spin-2 field in the fixed coordinate system 
such that the first axis (e x ) is aligned with <j) — 0, we obtain: 



±X'{n) 



I 



d 2 £ 

(27T) 2 



.X(l)e 



±2i<bp il*5 



(C5) 



where 4 + il y = Expanding ±X(1) = E{\) + iB(l), 

we can relate these rotationally invariant quantities to the 
Fourier transforms of the spin-2 field individual components. 



Cartesian map: 50 inversions 



10" 



o 



10" 



10 



-10 




Noise powerspectrum 



J I 



0.05 0.1 0.2 0.5 1.0 2. 

2n l/lmax 

Figure CI. The effect of noise on the reconstruction of the pow- 
erspectrum for a set of 50 realizations of the noise for the map 
1024 C/ in (^max = 1200). Note that at these scales, the spread in 
the recovered powerspectra for the different realizations is only 
visible above the cutoff frequency. 



In the case of weak lensing, we get the following flat sky 
limits: 



re(n) 



2 $(l)e l1 '" 



(C6) 



2 J (2tt) 2 

(7i±*ys)'(ft) « -\f (02^ 2 *(!) 

After identification, we thus get the limits for the operator 
J: 

J = (1,0) (C8) 
independently of the Fourier mode modulus. 

C2 SNR investigation in the plane 

Let us briefly investigate the effect of noise on the recov- 
ery of the k powerspectrum arising from the finite number 
of sources per unit area. For this purpose, let us consider 
the simplest setting corresponding to a cartesian map with- 
out mask which can therefore be inverted linearly following 
equations (122[1 - <I23I ). In this regime, the regularized solution 
is simply given in Fourier space by 



1 



flr 



+ 



^y) \^x ^y) 



, (C9) 



where g x , g y and k are the Fourier transform of the observed 
shear and convergence, and [i the penalty hyperparameter. 
In Figure [Cll we make use of the simulation 1024 Cf in , whose 
residues (after non linear inversion) are shown in Figure ([§]). 
Here 50 Monte Carlo realizations of the noise corresponding 
40 galaxies/Darcmin are averaged to produce an estimate 
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Expansion factor 

Figure Dl. the expected maximum uncertainty on particle posi- 
tions due to the method used to create the light cone as a function 
of the expansion factor. It is computed according to equation (ID 1 1 ) 
with a velocity v estimated to be 3 times the Virial velocity of 
the largest cluster in the simulation. 

of the corresponding errors. Clearly the shot noise remains 
small at all considered frequencies. cpp 



APPENDIX D: CONVERGENCE MAPS 

The inversion technique described in the main text was vali- 
dated using the mocks ext racted from the HORIZON-47T simu- 
lation (|Prunet et al.ll2008l ) . Let us briefly describe here how 
this simulation was used to generate mock slices and k maps. 

Dl Light cone generation 

The generation of a light cone during run time can be per- 
formed easily at each coarse time step of the simulation. 
Given a choice of the observer position in the simulation 
box, that we suppose here for simplicity to be at the ori- 
gin of coordinates, it is easy to select the particles that be- 
long to the slice in between redshifts Z2 < z\ corresponding 
to two successive coarse time steps: if (a;, y, z) are the co- 
moving coordinates of a particle, and d = \f x 2 + y 2 + z 2 
its comoving distance from the observer, we must have 
ddist(z2) < d < ddist(zi) for the particle to be selected, 
where ddist (z) is the comoving distance that a photon covers 
between redshift z and present time in the simulation box: 
rfdist = / cdt/a(t), where c is the speed of light and a the ex- 
pansion factor. The problem is that structures evolve during 
a coarse time step, so there are necessarily some discontinu- 
ities at the border between two successive light cone slices. 
These discontinuities are due to large scale motions of par- 
ticles plus their thermal velocity within dark matter halos. 
Given the large size of the simulation considered here, ther- 
mal motion within the largest cluster are expect to bring 
the most significant effects of discontinuity. For a particle 
with peculiar velocity v, the largest discontinuity to be ex- 
pected, i.e. the largest possible difference between expected 
and actual position of the particle is given by 

A = (i>/c)[ddiBt(*i) - d iist (z 2 )]. (Dl) 



In equation ()D1|) . we performed a linear Lagrangian approxi- 
mation, i.e. we neglected variations of the velocity of the par- 
ticle during the coarse time step. Using Press & Schechter 
formalism, or the improved formula of Sheth & Tormen 
(1999), the mass of the largest cluster in the Horizon simu- 
lation solves approximately the implicit equation 

Q. p c I?F[M ma3 ,{z),z\/M max (z) = 1, (D2) 

where p c is the critical density of the Universe and F is the 
fraction of mass in the Universe in objects of mass larger 
than M. Basically, this equation states that the mass in 
objects of mass larger than M is equal to M, which means 
that we are left with only one cluster of mass M, the largest 
detectable cluster in our cube of size L. We can compute 
F(M, z) with the usual formula, e.g. 

F(M, z) = f f(ji)dn, (D3) 

with v — 1.686/<t(M, z) where a(M,z) is the linear vari- 
ance at redshift z correspond ing to mass scale M, an d /(/i) 
is given by equation (10) of ISheth fc Tormenl (|l999T ). Per- 
forming these calculations, we find that the largest cluster 
at present time in a cube of size L = 2000/i _1 Mpc should 
have a typical mass of M max (z = 0) ~ 1.47 X 10 15 Af Q . With 
a standard Friend-of-friend algorithm using a linking pa- 
rameter 6 = 0.2, we find that the most massive halo de- 
tected in the simulation presents a somewhat larger mass, 
M = 5.4 x 1O 15 M0. Yet, in that rare events regime, we 
cannot expect our theoretical estimate to be more accurate. 
What matters, though, is the thermal velocity rather than 
the mass. Applying the Virial theorem, we have (e.g., Pea- 
cock, 1999) v 2 ~ GM max /R vir , with 

^nRl irPviT = Af max , p vir ~ 178fi /9c(l + z) 3 /Q(z) ' 7 , 

where Q(z) is the density parameter as a function of red- 
shift (fi(0) = fio)- These expressions are given in phys- 
ical coordinates hence the factor (1 + z) 3 in the expres- 
sion of pvir- This reads, at z — 0, v ~ 1570 km/s for 
M max (z = 0) ~ 1.47 x 1O 15 M . In the largest cluster of 
the simulation, the overall velocity dispersion is of the order 
of 2100 km/s, a slightly larger value that reflects the ac- 
tual value of the mass. To be conservative, we estimate the 
expected errors in equation <|D 1 1) with the Virial velocity 
rescaled by a factor 2100/1570, and with a further multipli- 
cation by a factor 3 to be in the 3<r regime. The correspond- 
ing maximal expected discontinuity displacement is shown 
in Mpc as a function of the expansion factor on Figure iDll 
As expected from the dynamically self-consistent calculation 
of the coarse time step (which is basically determined by a 
Courant condition using the velocity field), the comoving er- 
ror does not change significantly with redshift and remains 
below the very conservative limit of 200 kpc. Obviously, we 
expect in practice the errors brought by discontinuities to be 
in general much smaller than that, as for z = the present 
errors corresponds to unrealistic velocities as large as about 
6000 km/s! 

D2 From slices to k maps 

In the main text, the expression for k as a function of the 
density contrast in the simulation is given in equation (|28|l 
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in the geometric optic approximation. Let us rearrange this 
formula in a form that is more suited to integration over 
redshift slices in a simulation. 



2 ^ c JAz,, H &(Z) Ho 



where 
W 6 = 



dz r>(z)©(2,« 3 ) i \ / r dz 

Uh E{z) V&) ^[z)J 1 \L Zi W) 



is a slice-related weight, and the integral over the density 
contrast, 5, reads 

/">£/- \ V(simu) /iVpart(fpix,26) , 

Jax,, iVpart(smiu) V S p ^{z b ) 



where 

47T C 2 

is the comoving surface of the spherical pixel. Putting all 
together, we get the following formula for the convergence 
map: 

iVpix ( H \ 3 V(simu) ^ N pa , It (O p i x , z b ) 



2 4tt V c / A r P art(simu) ^ D 2 ( 2() 

Once the k map is available it is straightforward to build 
the corresponding g using equation (|17|l . 
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